QCD Dirac operator at nonzero chemical potential: lattice data and matrix model 
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Recently, a non-Hermitian chiral random matrix model was proposed to describe the eigenvalues 
of the QCD Dirac operator at nonzero chemical potential. This matrix model can be constructed 
from QCD by mapping it to an equivalent matrix model which has the same symmetries as QCD 
with chemical potential. Its microscopic spectral correlations are conjectured to be identical to 
those of the QCD Dirac operator. We investigate this conjecture by comparing large ensembles of 
Dirac eigenvalues in quenched SU(3) lattice QCD at nonzero chemical potential to the analytical 
predictions of the matrix model. Excellent agreement is found in the two regimes of weak and strong 
non-Hermiticity, for several different lattice volumes. 

PACS numbers: 12.38.Gc, 02.10.Yn 
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There has been a lot of recent interest in physical sys- 
tems described by non-Hermitian operators. Such oper- 
ators play a role in many areas of physics, e.g., S-matrix 
theory jlj] , dissipative quantum maps 2] , neural network 
dynamics 13], disordered systems with imaginary vector 
potential jj], and quantum chromodynamics (QCD) at 
nonzero density 0. In the present work, we are mainly 
interested in the last of these applications, but we expect 
the matrix model to be described below to be applicable 
to non-Hermitian operators in other physical systems as 
well, provided they are in the same symmetry class. 

QCD at nonzero density is important in a variety of 
physical situations, such as relativistic heavy- ion colli- 
sions or neutron stars. Considerable progress has been 
made in the last few years on the analytical side. For ex- 
ample, the regime of asymptotically large density is well 
understood |6|, and qualitative predictions for the QCD 
phase diagram could be derived on the basis of symme- 
try considerations However, quantitative results at 
physically relevant densities are still lacking at present. 
Unfortunately, lattice simulations of full QCD at nonzero 
chemical potential fj, are extremely difficult: the weight 
function is complex, and the numerical effort increases 
exponentially with the volume. A number of interesting 
new ideas have recently been investigated on the lattice 
side, e.g., reweighting along the critical Hne [31, combined 
expansions of weight function and observable [9j , analytic 
continuation from imaginary /i jlC|, and a factorization 
method for distribution functions of observables 0| . It is 
questionable, however, whether the techniques using real 
fi will allow us to approach the thermodynamic limit. 

Clearly, a better theoretical understanding of QCD at 
nonzero density is desirable. The Dirac operator is one 
of the central objects in QCD. Many observables can be 
expressed in terms of its eigenvalues and eigenvectors. 



While much is known about this eigenvalue spectrum at 
/i = (sec Rcf. 12] for a review), the situation at /i 7^ is 
less satisfying. The goal of the present work is to improve 
our understanding of the latter case. We concentrate on 
a particular matrix model for the QCD Dirac operator at 
nonzero fi and show that its analytical predictions for the 
distribution of small Dirac eigenvalues are in agreement 
with data from lattice gauge simulations. This statement 
holds in the two different regimes of weak and strong 
non-Hermiticity, to be defined below. The implications 
of these results are discussed in the conclusions. 

We start by presenting the matrix model and its pre- 
dictions [3] • The model constitutes a complex extension 
of the chiral Gaussian Unitary Ensemble (CUE) 0]. In 
terms of the complex eigenvalues Zj (j = 1, . . . , N), its 
partition function reads 
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for N f flavors of masses m/ (/ = !,..., iV/) in the sector 
of topological charge v. The parameter t€ JO, 1] measures 
the degree of non-Hermiticity and is related to by 



2 1 2 
/i = 1 — r 



(2) 



In the limit t — )■ 1 (or /i ^ 0), the eigenvalues are real 
and we are back to the chiral CUE. For t ^ the non- 
Hermiticity is maximal, and the model becomes a chiral 
extension of the Ginibre Ensemble . The relation (j^J 
follows from comparing the current model to the matrix 
model of Ref. at small /i. That model has the same 
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global symmetries as QCD and is defined by 
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The current model is equivalent to the model Q at 
the level of the partition function for small values of fi 
[l6| . (However, unlike the model lO, it is always in the 
phase with broken chiral symmetry.) So far it has only 
been possible to compute the microscopic spectral corre- 
lations (i.e., the correlations of the smallest eigenvalues 
on the scale of the mean level spacing) for model 
[l3| . and not for model It is conjectured that the 
two models, as well as QCD, are in the same universality 
class in the sense that they yield identical results for the 
microscopic spectral correlations. The analytical predic- 
tions for model can be derived at large N using the 
technique of orthogonal polynomials in the complex plane 
[l3|. All correlation functions have been obtained either 
for Nf = 0, or for Nf ^ "phase-quenched" massless 
flavors, replacing z^^^ \zj\'^^f . 

Two different large- TV limits have to be distinguished. 
In the limit of weak non-Hermiticity 0| , the product 
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is kept fixed. This corresponds to taking the volume 
V (x N to infinity such that Vfj,^ is fixed. The result of 
the current model for the density p{z) = (Iljli ^{^"^j)) 
of small Dirac eigenvalues in this limit reads at Nf = 
(corresponding to the quenched limit in our lattice data) 
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where J denotes the Bessel function and the eigenvalues 
have been rescaled according to ^ = V2Nz, resulting in 
the customary level spacing of tt. The limit of strong 
non-Hermiticity, with iV — > 00 at fixed re [0, 1), leads to 



Pstrong(e) = V2^|?|exp(-|e|')/|.|(|C|' 



(6) 



where / denotes the modified Bessel function. The rescal- 
ing ^ = N/{l — T^) z results in a level spacing indepen- 
dent of r. We stress that the existence of these two dif- 
ferent scaling regimes is a prediction for the lattice, and 
we will identify these two regimes in the data below. 

We now turn to the details of the lattice simulations 
and discuss some of the concerns arising from our choices 
of operator and simulation parameters. We use the stag- 
gered Dirac operator, given at 7^ in terms of SU(3) 



jauge fields U and staggered phases r? by 

i^—x,y.z 



(7) 



The lattice spacing has been set to unity. We denote 
its eigenvalues by iXk with Afc real (complex) for /i = 
(/i 7^ 0). The reason to prefer the staggered formulation 
is that (a) the Wilson operator breaks chiral symmetry 
explicitly and has complex eigenvalues even at /i = 
and (b) Ginsparg- Wilson-type operators are much more 
expensive to compute. This is a serious issue here because 
we need many configurations (see Table HJ. 

The gauge field configurations were generated in the 
quenched approximation, corresponding to Nf = 0. At 
/i ^ the Nf — > limit is subtle. In the current model, 
this limit can be taken in different ways: (i) use Eq. 
and let Nf ^ at the end of the calculation, (ii) use 
Eq. "phase-quenched", i.e. with |zjp^-'' in the weight, 
and let Nf at the end of the calculation, or (iii) set 
Nf = in Eq. at the beginning of the calculation, (ii) 
and (iii) yield identical results, given in Eqs. H5I6|I . and 
these results agree with quenched lattice data at /i ^ 0, 
as shown below. This is consistent with the fact that 
quenched QCD at /i 7^ is the Nf limit of a theory in 
which (det D)^f is replaced by | det D\^f 5J. Taking the 
Nf — > limit according to (i) will lead to a different (yet 
unknown) result which would not describe the quenched 
lattice data. At Nf > 0, however, we expect Eq. to be 
valid for unquenchcd QCD, since it contains the correct 
phase of the fcrmion determinant. This expectation is 
strengthened by the fact that matrix models of this type 
have already been successful in describing unquenched 
QCD at /I 7^ 0,0|. (We also note that matrix models 
describe unquenched lattice data at = [13 ■) 

In the simulations we used /3 = — 5.0, which cor- 
responds to the strong-coupling regime and is far from 
the continuum limit. Nevertheless, working at such a 
low value of (3 is both convenient and legitimate for the 
present purpose, for the following reason. The micro- 
scopic spectral correlations are described by a random 
matrix model only below the so-called Thouless energy, 
which is the boundary of the regime in which the zero- 

TABLE I: Summary of simulation parameters (/? — 5.0). 



V 


^J■ 


level spacing d 


no. of confi 


6* 


0.006 


1.98(2) • 10"^ 


17,000 


6* 


0.03 


1.55(3) • 10"^ 


20,000 


6* 


0.2 


6.83(6) • 10"^ 


20,000 


8* 


0.003375 


6.30(3) • 10"* 


20,000 


8* 


0.2 


3.85(3) • 10"^ 


20,000 


10* 


0.00216 


2.57(4) • 10"* 


4,000 


10* 


0.2 


2.46(4) ■ 10"^ 


4,000 
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momentum modes dominate the partition function of the 
low-energy effective theory in a finite volume [l^ . The 
Thouless energy is a function of both /3 and the lattice 
volume V. If jS is increased at fixed V, fewer Dirac eigen- 
values are described by the matrix model. Increasing V 
works in the opposite direction. Thus, for small /3 we 
can test the matrix model on relatively small lattices. At 
larger values of /3 we simply need to increase the lattice 
volume, which, however, is inconvenient numerically. 

At this small value of /?, the staggered Dirac opera- 
tor does not have exact zero modes even if the underly- 
ing gauge field has nonzero topological charge, because 
the would-be zero modes are shifted by an amount pro- 
portional to the square of the lattice spacing . This 
amount is much larger than the mean level spacing near 
zero, and thus the would-be zero modes are completely 
mixed with the nonzero modes. We account for this by 
setting J/ = in Eqs. H1I5I6() . (This "disease" of staggered 
fermions can be overcome by going to very small lattice 
spacing or by using the overlap operator H^.) 

Our simulation parameters are summarized in Table 
In the regime of strong non-Hermiticity, we used a con- 
stant value of /i = 0.2, whereas for weak non-Hermiticity, 
we varied such that the product n'^V is fixed. The 
gauge fields were generated using a combined overre- 
laxation and Metropolis algorithm, written by P.E.L. 
Rakow. On the B'* lattice, the complete eigenvalue spec- 
trum can be c omp uted on a single PC rather quickly us- 
ing LAPACK |23j. For the larger lattices we switched to 
ARPACK ji^l and computed only the 100 eigenvalues of 
smallest absolute magnitude with positive real part (the 
eigenvalues come in pairs ±Afc). Since ARPACK is very 
fast at computing the largest eigenvalues, we inverted the 
Dirac operator prior to feeding it to ARPACK, using the 
sparse LU solver UMFPACK 2^. 

Because we are only interested in the small eigenvalues 
of the Dirac operator, global unfolding of the spectrum 
|2(tI | is not necessary; we simply rescale the eigenvalues 
by a constant determined from the mean level spacing d 
of the data near zero. Note that d oc 1/V in the weak 
limit and d oc 1/VV in the strong limit, respectively. 

We first present our results for weak non-Hermiticity. 
Since three-dimensional plots are hard to read, we instead 
show cuts along the real and imaginary axes. In Fig. ^ 
the data for lattice size 6^ and fi = 0.006 are plotted ver- 
sus Eq. IS)). There is no free fit parameter; the data have 
been rescaled according to ^ = irX/d, with d ocl/V given 
in Table (note that d depends on the lattice spacing, 
i.e. it would change with /?). At weak non-Hermiticity 
d can be obtained in the same way as for real eigenval- 
ues: Because of the smallness of their imaginary part 
the eigenvalues can still be ordered with respect to their 
real part so that the level spacing is defined unambigu- 
ously. The very same level spacing d is used to determine 
^2 = {tt/V2) ij?/d from Eq. (gj) and dmodci = tt/{'\/2N), 
leading to a = 0.20 for use in Eq. 10 • Within our sta- 
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FIG. 1: Density of small Dirac eigenvalues for y = 6* and 
= 0.006, cut along the real axis (left) and parallel to the 
imaginary axis at the first maximum (right). The histogram 
represents lattice data, and the solid curve is Eq. @. 



tistical accuracy, we obtain excellent agreement of lattice 
data and analytical prediction without any free fit param- 
eter. In Fig. 121 we repeat the same analysis for lattice size 
8^ and /i = 0.003375, chosen to keep V ji^ — c? constant 
in order to test the scahng predicted by Eq. Here 
we have used the same value of a = 0.20 and again find 
excellent agreement. This value agrees within 0.3% with 
the value of a determined independently by rescaling /i^ 
with the level spacing d from the S'' data. Similar results 
are obtained from the analysis of the 10'' data which we 
do not display here because of limited statistics. 

Next we turn to the limit of strong non-Hermiticity. 
Here, the eigenvalues are no longer localized close to the 
real axis and spread into a two-dimensional domain in the 
complex plane. Therefore we have to modify the deter- 
mination of the mean level spacing d. Nearest neighbors 
are now defined as having the smallest geometric distance 
between them. As can be seen from Table the level 
spacing d oc 1/VV is very different from the weak limit. 
The signature of the chiral ensemble compared to the 
Ginibrc ensemble Is^ is a "hole" at the origin resulting 
from the level repulsion. (The latter ensemble applies 
only in the bulk of the spectrum 26].) To compare the 
lattice data with Eq. lO, we rescale the eigenvalues ac- 
cording to ^ = cA/d, where c is a constant related to the 
level spacing of our model in the strong non-Hermiticity 
limit. We currently do not have a theoretical result for 
c and therefore determined c — 0.82(5) by a fit to the 
data 113. In Figs. 01 and H we plot the data for ^ = 0.2 
and V — 6'^, 8^ versus the prediction of Eq. ©. Good 
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FIG. 3: Density of small Dirac eigenvalues for V — 6'^ and 
fi — 0.2, cut along the real (left) and imaginary (right) axes. 
The histogram represents lattice data, and the solid curve is 
the prediction of Eq. 

agreement is found along both the real and imaginary 
axis (also for V = 10'*, not shown), confirming the rota- 
tional invariance of the microscopic density. Note, how- 
ever, that there is an important difference between the 
cuts in the real and imaginary directions. While the mi- 
croscopic density, Eq. 10, is rotationally invariant, the 
data spread macroscopically into a thin ellipse. Along 
the imaginary axis the spectrum ends at « ±3.8 in our 
units, while along the real axis it extends up to « ±270. 
For that reason the last part of the histograms along the 
imaginary axis in Figs. |3| and ^ may no longer be in the 
microscopic regime. In the data for 1/ = 6^ at interme- 
diate = 0.03, the microscopic and macroscopic scales 
no longer separate clearly, and therefore neither of the 
Eqs. (O and lO apply. 

In conclusion, we have identified two different regimes 
in the behavior of complex eigenvalues of the QCD Dirac 
operator in the domain where chiral symmetry is broken. 
Our lattice data confirm the predictions of random ma- 
trix theory quantitatively, both at weak and at strong 
non-Hermiticity. Matrix models thus provide a detailed 
theoretical understanding of the properties of complex 
Dirac eigenvalues in a new regime with non-vanishing 
chemical potential, including the correct scaling of eigen- 
values and chemical potential with the volume. Our find- 
ings may have algorithmical implications, since it is typ- 
ically the low-lying Dirac eigenvalues which determine 
the numerical effort in lattice simulations. We wish to 
emphasize that, although our conclusions are based on 
quenched simulations using staggered fermions, the pre- 
dictions of the matrix model could also be tested in un- 
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FIG. 4: Same as Fig.glbut for = S"* and p = 0.2. 



quenched simulations and in sectors of nontrivial topo- 
logical charge. This will be the subject of future work. 

This work was supported by the DFG (G.A.) and in 
part by DOE grant No. DE-FG02-91ER40608 (T.W.). 

Note added in proof. — The microscopic density of the 
model Eq. Q was recently computed in the weak non- 
Hermiticity limit [23. The result deviates very slightly 
from our Eq. (O, but this difference cannot be resolved 
by our data. 
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Erratum: QCD Dirac Operator at Nonzero Chemical Potential: 
Lattice Data and Matrix Model 
[Phys. Rev. Lett. 92, 102002 (2004)] 
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PACS numbers: 12.38.Gc, 02.10.Yn, 99.10. Cd 

In Figs. 1 and 2 we compared lattice data to the analytical prediction of Eq. (5), which pertains to the weak non- 
Hermiticity limit of the model in Eq. (1). For such a comparison, two independent parameters must be determined: 
the mean level spacing d and the dimcnsionless parameter a appearing in Eq. (5). These two parameters can be 
expressed in terms of two low-energy constants, the chiral condensate T, at fi — and the pion decay constant Fj^, 
by the relations d = n/J^V and = 2^^F^V 1]. Here, V is the physical volume. We erroneously determined as 
(/ua)^7r/-\/2fia, where a is the lattice spacing. This amounts to setting (FTra)^/I]a^ = 2^^/^. Instead, F-^ and thus a 
should be obtained either by an independent measurement of or by a fit of the data to Eq. (5). Performing the 
latter for the data sets in Figs. 1 and 2 yields a value of a = 0.200(2), which is consistent with the value of a = 0.20 
used in the Letter. Thus, Figs. 1 and 2 as well as the conclusions of the Letter remain unchanged. 

The fitting of a to the lattice data constitutes a new method for obtaining F^^ through the relation quoted above, see 
also Ref. , but we refrain from applying this method here since the lattice data were obtained in the strong-coupling 
region. A related method to determine using isospin chemical potential was introduced recently in Ref. 0]. 

We thank J.C. Osborn for pointing out the error described above and P.H. Damgaard for useful discussions. 
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